%This code produces 6 plots that present the comparative statics wrt sigma, and
%baseline plots. 

%Last update 03/01/2023

clear all
close all
clc


myopt = optimset('Display','Off','MaxFunEvals',1e1000000,'TolFun',1e-20,'TolX',1e-30);

warning off
datadir = '..\..\Data\Final\Figures8_9_10\';
figdir = '..\..\Plots\';
load([datadir 'sigma_data_S100J50_alpha_0001_ADJalpha']);


%%
vec_q=[grid_q(1) 1 grid_q(end)];
q3=vec_q(3);
q2=vec_q(2);
q1=vec_q(1);

for qq=1:length(vec_q)
    q=vec_q(qq);
    for m=1:M
        for s=1:S
            n_s = pchip(grid_q,n_mat_save(:,s,m));
            sales_s = pchip(grid_q,sales_mat_save(:,s,m));
            plants(qq,s,m) = ppval(n_s,q);
            sales(qq,s,m) = ppval(sales_s,q);
        end
    end
end

fig1 = figure(1);
qh=plot(log(grid_income),plants(3,:,1),'linewidth',2.5,'Color',[0 0 0]);
hold on
qh=plot(log(grid_income),plants(2,:,1),'linewidth',2.5,'Color',[0 0 1]);
hold on
qh=plot(log(grid_income),plants(1,:,1),'linewidth',2.5,'Color',[91, 207, 244] / 255);
th=legend("$q=$"+q3+"","$q=$"+q2+"","$q=$"+q1+"",'Location','northwest');
set(th,'Fontsize',9,'Interpreter','Latex')
th=ylabel("Measure of Plants of Firm $j$ in $s$, $n_{js}$");
set(th,'Fontsize',12,'Interpreter','Latex')
th=xlabel("Log of Income in $s$, $\ln I_s$");
set(th,'Fontsize',12,'Interpreter','Latex')
axis tight
saveas(fig1,[figdir 'Figure8_left'],'png');
%%
fig2 = figure(2);
qh=plot(log(grid_income),plants(3,:,1),'linewidth',2.5,'Color',[0 0 0]);
hold on
qh=plot(log(grid_income),plants(2,:,1),'linewidth',2.5,'Color',[0 0 1]);
hold on
qh=plot(log(grid_income),plants(1,:,1),'linewidth',2.5,'Color',[91, 207, 244] / 255);
hold on
qh=plot(log(grid_income),plants(3,:,M),'linewidth',2.5,'Color',[0 0 0],'linestyle',':');
hold on
qh=plot(log(grid_income),plants(2,:,M),'linewidth',2.5,'Color',[0 0 1],'linestyle',':');
hold on
qh=plot(log(grid_income),plants(1,:,M),'linewidth',2.5,'Color',[91, 207, 244] / 255,'linestyle',':');
th=legend("$q=$"+q3+", $\sigma=$"+run_vec(1)+"","$q=$"+q2+", $\sigma=$"+run_vec(1)+"","$q=$"+q1+", $\sigma=$"+run_vec(1)+"","$q=$"+q3+", $\sigma=$"+run_vec(M)+"","$q=$"+q2+", $\sigma=$"+run_vec(M)+"","$q=$"+q1+", $\sigma=$"+run_vec(M)+"",'Location','best');
set(th,'Fontsize',9,'Interpreter','Latex')
th=ylabel("Measure of Plants of Firm $j$ in $s$, $n_{js}$");
set(th,'Fontsize',12,'Interpreter','Latex')
th=xlabel("Log of Income in $s$, $\ln I_s$");
set(th,'Fontsize',12,'Interpreter','Latex')
axis tight
saveas(fig2,[figdir 'Figure9_left'],'png');
%%
fig3 = figure(3);
for m=1:M
    for j=1:length(vec_q)
        for s=1:S
            if sales(j,s,m)<=0
                sales_plot(j,s,m)=0.00000001;
            else
                sales_plot(j,s,m)=sales(j,s,m);
            end
        end
    end
end

qh=plot(log(grid_income),log(sales_plot(3,:,1)),'linewidth',2.5,'Color',[0 0 0]);
hold on
qh=plot(log(grid_income),log(sales_plot(2,:,1)),'linewidth',2.5,'Color',[0 0 1]);
hold on
qh=plot(log(grid_income),log(sales_plot(1,:,1)),'linewidth',2.5,'Color',[91, 207, 244] / 255);
th=legend("$q=$"+q3+"","$q=$"+q2+"","$q=$"+q1+"",'Location','northwest');
set(th,'Fontsize',9,'Interpreter','Latex');
th=ylabel("Log of Sales in $s$");
set(th,'Fontsize',12,'Interpreter','Latex');
th=xlabel("Log of Income in $s$, $\ln I_s$");
set(th,'Fontsize',12,'Interpreter','Latex');
axis([log(grid_income(1)) log(grid_income(end)) -2 max(log(sales_plot(end,:,end)))])
saveas(fig3,[figdir 'Figure8_right'],'png');
%%
fig4 = figure(4);
qh=plot(log(grid_income),log(sales_plot(3,:,1)),'linewidth',2.5,'Color',[0 0 0]);
hold on
qh=plot(log(grid_income),log(sales_plot(2,:,1)),'linewidth',2.5,'Color',[0 0 1]);
hold on
qh=plot(log(grid_income),log(sales_plot(1,:,1)),'linewidth',2.5,'Color',[91, 207, 244] / 255);
hold on
qh=plot(log(grid_income),log(sales_plot(3,:,end)),'linewidth',2.5,'Color',[0 0 0],'linestyle',':');
hold on
qh=plot(log(grid_income),log(sales_plot(2,:,end)),'linewidth',2.5,'Color',[0 0 1],'linestyle',':');
hold on
qh=plot(log(grid_income),log(sales_plot(1,:,end)),'linewidth',2.5,'Color',[91, 207, 244] / 255,'linestyle',':');
th=legend("$q=$"+q3+", $\sigma=$"+run_vec(1)+"","$q=$"+q2+", $\sigma=$"+run_vec(1)+"","$q=$"+q1+", $\sigma=$"+run_vec(1)+"","$q=$"+q3+", $\sigma=$"+run_vec(M)+"","$q=$"+q2+", $\sigma=$"+run_vec(M)+"","$q=$"+q1+", $\sigma=$"+run_vec(M)+"",'Location','northwest');
set(th,'Fontsize',9,'Interpreter','Latex');
th=ylabel("Log of Sales in $s$");
set(th,'Fontsize',12,'Interpreter','Latex')
th=xlabel("Log of Income in $s$, $\ln I_s$");
set(th,'Fontsize',12,'Interpreter','Latex')
axis([log(grid_income(1)) log(grid_income(end)) -2 max(log(sales_plot(end,:,end)))])
saveas(fig4,[figdir 'Figure9_right'],'png');
%% Markets
t = @(x) (exp(sqrt(1/phi)*x)).^(1-eps) ;
kappa_prime_zero = integral(@(u)(t(u).*u*2*pi),0,Inf);
R_function = @(income) (exp(rent_1*(log(income)).^rent_2));

income_l=zeros(J,M);
income_h=zeros(J,M);
for m=1:M
    sigma=run_vec(m);
    Zs_spline = spline(grid_income,Zs_vec_save(m,:));
    Zs_function = @(income) ppval(Zs_spline,income);
    z_function = @(q,N,sigma) q.^(eps-1).*(exp(-(1/sigma).*N)).^(eps-1) ;
    for j=1:J
        if sum(n_mat_save(j,:,m)>0)>0
            q=grid_q(j);
            lambda=lambda_J_save(j,m);
            N=N_J_save(j,m);
            zj=z_function(q,N,sigma);
            [value,ubics]=find(n_mat_save(j,:,m)>0);
            income_l_aux=grid_income(ubics(1));
            income_h_aux=grid_income(ubics(end));
            solve_fun = @(inc) (((inc/(eps-1))/Zs_function(inc))*zj*kappa_prime_zero-R_function(inc)+lambda)^2;
            if income_l_aux==1
                income_l(j,m)=1;
            else
                income_l_aux2=grid_income(ubics(1)-1);
                income_l(j,m) = fmincon(@(n) solve_fun(n),income_l_aux/2+1/2,[],[],[],[],income_l_aux2,income_l_aux,[],myopt);
            end
            if ubics(end)==S
                income_h(j,m)=grid_income(end);
            else
                income_h_aux2=grid_income(ubics(end)+1);
                income_h(j,m) = fmincon(@(n) solve_fun(n),income_h_aux+20,[],[],[],[],income_h_aux,income_h_aux2,[],myopt);
            end
        else
            income_l(j,m)=1;
            income_h(j,m)=1;
        end
        clear value ubics
    end
end

fig5 = figure(5);
qh=plot(log(grid_income),log(grid_income/(eps-1)./Zs_vec_save(1,:)),'linewidth',2.5,'Color',[0 0 0]);
hold on
qh=plot(log(grid_income),log(grid_income/(eps-1)./Zs_vec_save(end,:)),'linewidth',2.5,'Color',[0 0 0],'linestyle',':');
th=legend("$\sigma=$"+run_vec(1)+"","$\sigma=$"+run_vec(end)+"",'Location','northwest');
set(th,'Fontsize',9,'Interpreter','Latex')
th=ylabel("Log of Profitability of Location $s$, $\ln x_s$");
set(th,'Fontsize',12,'Interpreter','Latex')
th=xlabel("Log of Income in $s$, $\ln I_s$");
set(th,'Fontsize',12,'Interpreter','Latex')
axis tight
saveas(fig5,[figdir 'Figure10_right'],'png');
%%
fig6 = figure(6);
qh=plot(log(grid_q),lambda_J_save(:,1),'linewidth',2.5,'Color',[0 0 0]);
hold on
qh=plot(log(grid_q),lambda_J_save(:,end),'linewidth',2.5,'Color',[0 0 0],'linestyle',':');
th=legend("$\sigma=$"+run_vec(1)+"","$\sigma=$"+run_vec(end)+"",'Location','northwest');
set(th,'Fontsize',9,'Interpreter','Latex')
th=ylabel("Marginal Cost of an Extra Plant for Firm $j$, $\lambda_j$");
set(th,'Fontsize',12,'Interpreter','Latex')
th=xlabel("Log of Firm $j$ Productivity, $\ln q_j$");
set(th,'Fontsize',12,'Interpreter','Latex')
axis tight
saveas(fig6,[figdir 'Figure10_left'],'png');
